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Abstract 

We present solution of self-consistent equations for the N 3 LO nuclear energy 
density functional. We derive general expressions for the mean fields ex- 
pressed as differential operators depending on densities and for the densities 
expressed in terms of derivatives of wave functions. These expressions are 
then specified to the case of spherical symmetry. We also present the com- 
puter program hosphe (vl.OO), which solves the self-consistent equations by 
^C) ', using the expansion of single-particle wave functions on the spherical har- 

monic oscillator basis. 
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Nature of problem: 

The nuclear mean-field methods constitute principal tools of a description of nu- 
clear states in heavy nuclei. Within the Local Density Approximation with gra- 
dient corrections up to N 3 LO [1], the nuclear mean-field is local and contains 
derivative operators up to sixth order. The locality allows for an effective and fast 
solution of the self-consistent equations. 
Solution method: 

The program uses the spherical harmonic oscillator basis to expand single-particle 
wave functions of neutrons and protons for the nuclear state being described by 
the N 3 LO nuclear energy density functional [1]. The expansion coefficients are 
determined by the iterative diagonalization of the mean-field Hamiltonian, which 
depends non-linearly on the local neutron and proton densities. 
Restrictions: 

Solutions are limited to spherical symmetry. The expansion on the harmonic- 
oscillator basis does not allow for a precise description of asymptotic properties of 
wave functions. 
Running time: 

50 sec. of CPU time for the ground-state of 208 Pb described by using the Nq = 50 
maximum harmonic-oscillator shell included in the basis. 
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[1] B.G. Carlsson, J. Dobaczewski, and M. Kortelainen, Phys. Rev. C 78, 
044326 (2008). 
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LONG WRITE-UP 



1. Introduction 

The nuclear mean-field methods constitute principal tools of a descrip- 
tion of nuclear states in heavy nuclei [1(1. Their applicability to nuclei is 
interpreted within the Hohenberg-Kohn [2| and Kohn-Sham theorems in- 
volving the nuclear energy density functional (EDF). Recently, we formulated 
the N 3 LO nuclear EDF with gradient corrections up to sixth order jij. The 
present study presents practical formulation of the method, allowing for a 
solution of the corresponding self-consistent equations. We also present the 
computer program hosphe (vl.OO), which solves the self-consistent equations 
by using the expansion of single-particle wave functions on the spherical har- 
monic oscillator (HO) basis. 

The paper is organized as follows. In Section (2], we present concise review 
of the method. In Section (3], we give general forms of the N 3 LO potentials, 
fields, and densities, which are then in Section H] specified to the case of 
spherical symmetry. Sections loT-fTUl describe the structure, installation, and 
test runs of the code hosphe (vl.OO), and Section [TT1 concludes our study. 

2. Overview of the method 

This introductory section is intended as a guide to the subsequent sec- 
tions, where more detailed derivations and results are presented. Here, we 
use abbreviated notations so as to give a brief outline of the method, while 
referring the reader to the following sections for details. 

The mean-field eigenvalue equation is obtained by considering the vari- 
ation of the energy with the condition of the single-particle wavefunctions 
(pi (r, a) to be normalized to unity, 



5^ " a 



(i) 



The potential energy is expressed as the EDF of Ref. [4|, that is, 
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where the grouped indices, such as the greek indices a = {n a L a v a J a } or 
P = {npLpVpJp) and the roman indices a = {m a I a }, denote all the quantum 
numbers of the local primary pp(r) and secondary p a ,a,j{r) densities (!]. In 
Eq. (j2J), T^ a {r) denotes terms of the functional, 



T a,a = [P(3pa,a] = 



Pd [Dapa\.j B 



Png 



(3) 



Cf a denote the coupling constants, D a denote the higher order derivative 
operators [H, and the sum runs over all terms of the functional. Although 
not shown explicitly, the sum contains both isoscalar and isovector terms. 
At present we have neglected neutron-proton mixing, which means that only 
the T z — component of the isovector densities are present. The convention 
adopted is such that the isovector contributions to the densities are taken as 
neutron minus proton densities and isoscalar densities are sums of neutron 
and proton densities. 

For each term of the functional, the variation with respect to the local 
densities, followed by the integration by parts and recoupling gives, 



5 / TLd 3 r 



[D a Pa)j a 



J o 



l8pa[D a p p ) Ja ] )d*r, (4) 



cf. Eqs. ( |43|) and (jHJ). The primary densities can be expressed in terms of 
the single-particle wavefunctions as: 



PnLvJ 



(5) 



The higher order derivative operators K n i jij are built by coupling the 
relative- momentum operators k = ^ (V^ — V^*), where we have used the 
subscripts to indicate on which function the operators act upon. This allows 
us to perform the variation with respect to 4>*, 



J a. 



[D aP , 



+ (-1) 



J Q — Jot 



[K a L a Vv a <j)i\ j0 [DaP/3] 



(6) 



After performing the variation, the integral above was partially integrated 
so that the derivatives would not act on the variation of <fi*. Therefore, 
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the K' nL operators are built by coupling the relative-momentum operators 
k' = Yi (V + , where V acts on all the functions of position standing to 
the right. 

The operator on the right-hand-side of Eq. ([6]) is a formal expression 
for the mean-field operator. All what remains to be done is to disentangle 
the gradients and V from one another - this procedure is performed 
in Eqs. below. Finally, the mean-field operator h (p) acquires the 
form: 

(7) 



TJ 



where the differential operators -D„ 7 l 7 and Pauli matrices a Vi act on the 
single-particle wave functions, and the potentials U 7 (r) are linear combina- 
tions of the secondary densities: 



^7 — ^ ^a,aXa,t^ [Ai 
aa(3;d8 



(8) 



The coefficients Xa'a^ can ^ e derived by using the recoupling rules pre- 



sented in Section 13.31 An alternative method, which was also used when 
building the code hosphe (vl.OO), was to construct the fields by starting 
from Eq. (jBJ) and putting them equal to those of Eq. ([7j) . This gives a linear 
system of equations that can be solved for the unknown coefficients xfaV 
At N 3 LO, only 1494 such coefficients are needed, so they can easily be pre- 
calculated and stored. 

It is now clear, that the key operators in the mean field are given by 



(9) 



and their matrix elements in the single-particle basis <pi (r, a) read, 



d<5,7 



•VyI&') 



J2<i>Ur,a) \p^} [D 



r'Ov r, a 



d 3 r. (10) 



Then, the mean-field matrix elements can be written as the following sum: 
(</H\h(p)\<k) = £ Cij&Mr ( n ) 

aa/3;d<57 
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Matrix elements in a spherical basis are derived in Sections 14.41 and 14.51 

When constructing potentials flS]), we need expressions to calculate all 
secondary densities. These can be written as [see Eqs. (jT6"]) -(l 



Pd,5,JM = [Daps] 



JM 



E*bb',W bb',W 
n d,6,J Pv s JMi 



bb'W 



with 



bb',W I n 
Pv s JM 1*% 



n(i)n( 2 ) 



Pv s [ri,r 2 ) 



(12) 



(13) 



JM ) r=r2=vi 

where the superscripts on the derivative operators indicate on which coordi- 
nate they act. The coefficients A can be obtained by explicitly constructing 
the left- and right-hand sides of Eq. (|T2"|) , which gives a linear system of equa- 
tions in derivatives of the density matrix that can be solved for the unknown 



coefficients A 



bb',W 
d,S,J 



At N 3 LO, only 3138 such coefficients are needed, so they 



can easily be precalculated and stored. In Section 13.61 we also show how to 
derive these coefficients by using the recoupling rules and in Section 14.21 we 
give the expressions for densities in the spherical HO basis. 



3. General forms of the N 3 LO potentials, fields, and densities 

3.1. Building blocks 

We begin by recalling definitions that are used to construct operators and 
densities in the spin and position space. The basic building blocks are given 
as in Eqs. (8)-(9) of Q, i.e., 

O~v=0,0 
C«=1,^={-1,0,1} 
Vi^={-l,o,l} 
fcl,^={-l,0,l} 

where k is the relative momentum operator: 

fc = ^(V-V). (18) 

All possible N 3 LO differential operators -D n LM, which can be built of gradi- 
ents ( Fl6|) . are given in the Table I of Ref. [4|, where n is the order of the 



^ (14) 

~ l { 72 ~ i(J y} ' a ^ 73 ( a * + i(J y} } > ( 15 ) 

-i {73 (V, - zV„) , V„ f 2 (V, + iV y )} (16) 

—% *f^g [k x — iky) , k z , (k x + ifcj,)! , (17) 
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= QA 




for 


Q CT = +1, 


n. 


= <2v< 




for 


Qv = -l, 




= g fc ( 




for 


Q k = +1. 



operator and L is its rank with magnetic projection M. Exactly in the same 
way, in Ref. [4J we defined the operators K nLM , which are spherical tensors 
built of the relative momentum operators k (fTTj) . 

Hermitian-conjugation properties of the building blocks read: 

(19) 
(20) 
(21) 

For any pair of commuting operators Ax^ and that have the following 
hermitian-conjugation properties: 

K = Qa(-1) X -»A x ^, (22) 
= QBi-lf-'Sw, (23) 

(24) 

the operator Clm built by the angular momentum coupling, 

C L m = [A\B X '}lm = 22 Cx^'p'AxuBx'n', (25) 

behaves under the hermitian conjugation as: 

Cl M = Qc(-1) L ~ M C l ^m for Qc = QaQb- (26) 
As a consequence, we have 

Klm= (-1)" {-lf- M D nL ,- M , (27) 
Klm = {-\) L ~ M K nU _ M . (28) 

We note that the gradient operators (TTB]) and (ITTj) obey the Biedenharn- 
Rose phase conventions of 

V;„ = (-l) 1 -^!,-^, (29) 
K,= - (-l) 1 "^!,-/*. (3°) 

^LM = (-l) i - M J DnX,-M, (31) 
^nLM = (-I)XlM, (32) 
^nLM = (-l) n D nLM i (33) 



which gives 
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and 

KLM = (-lT(-l) L - M KnL,-M, (34) 

Klm = (-ITKlm, (35) 
Klm = (-l) n K nLM , (36) 

where superscript T denotes the transposed operator. 
3.2. Potentials 

The potential energy to be varied over the wave functions is given in 
Ref. m and reads 



£ = J d 3 rH(r), (37) 

for 

n(r) = C ^nii'j[Pn'L^Ar)[D mlPnLvJ (r)l 



J' JO 



n'L'v'.l' 
ml ,nLv J 



^ CmI,nLvJ ( y/lj'+l Pn> L>v> J' M>(r)[D mI p nLvJ {r)}j> _ M t l?>%) 



V2J 

n'L'v'j'M' 
ml ,nLv J 



where C^f^vj are ^ e coupling constants and p n Lvj( r ) are the primary 
densities: 

pnLvj{r) = {[K nL p v (r, r')]j} rl=r , (39) 

which are built by acting with the relative momentum operators K n iM on 
the scalar (v = 0) and vector (v = 1) non-local densities: 

Pv„(r,r') = ^p(r(T,rV) (ff'K». (40) 

era' 

Note that the sum in Eq. (1381) runs over the indices ordered in a specific way, 
defined in Ref. [1, 0] , namely, 

{n'L'v'f} < {nLvJ}. (41) 

We first vary £ over the densities and then, in Section 13.31 we vary den- 
sities over the wave functions, that is, we begin with 

58= f d 3 r5H(r) = V / d 3 r- — — {r)5p n , Ur ,j, M ,. (42) 

J n'MJ'W J d Pn'L>v'J'M> 
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An explicit variation over the primary densities under the differential oper- 
ators D mI can be avoided by first integrating by parts and recoupling. The 
recoupling within a scalar [6( is simple, namely, 

[AABjCjUo = (-ly'+^iCjiBjAMo. (43) 

Hence, the integration by parts gives: 

n(r) = Yl (-^^"^^^WI^^w^WJjIo 

n'L'v' J' 
ml ,nLvJ 

= E i-^'^C^LvJ [PnLvj(r)[D mlPn , Uv ,j,(r))j}o 

n'L'v' J' 
ml ,nLvJ 

= J2 (- 1 ) J_J '^XLVj'[Pnavj'(^)[^/PnL,j(r)]^] , (44) 

n'L'v' J 1 
ml ,nLvJ 

where we have used Eq. fl3"3"|) . then we changed the names of indices, and we 
also used the fact that m + / is even. 

Therefore, the variation under the differential operators D m i only gives 
the transposition of indices {n'L'v' J'} «-> {nLvJ} and the phase. The com- 
plete variation of the energy then reads: 

5£ = E / d3r l S Pn'L'v'j'Un>L'v'j'(r)] , (45) 

n'L'v' J' J 

where we defined potentials 

U n 'L'v'j'M'(r) =E \^mi,ndj + ( -1 ) J ~ J C^wj'} [D m ip n Lvj{r)] jw(46) 

mI,nLv J 

Note that because of the ordering (I4TI) . in Eq. (f4*6|) either the first or the 
second term is non-zero (for {n'L'v' J'} 7^ {nLv J}), or both terms add up to 
2C$£& (for {n'LVJ'} = {nL^J}). 

5.5. FieZds 

Now we are in a position to perform the variation of densities over the 
wave functions. To this end, we assume that the non-local density matrices 
in Eq. (HUj) have the general form of 

p(ra,r'a') = M™)Mr V). (47) 

i 
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This form allows us to carry out the derivation for several important specific 
cases simultaneously. Namely, for the standard HF case one has: 

V*(rV) = #(rV), (48) 

and the sum runs over the occupied states only, i — 1, . . . , A. Similarly, for 
the BCS case, or for the HFB case in the canonical basis, one has: 

i(rV)=*V), (49) 

where vf are the occupation factors and the sum runs over the pairing win- 
dow. For transition densities pertaining to the symmetry restoration, one 
has: 

^(rV) = £<V(a)#(r V,«), (50) 

3 

Where ^(rV',a) = R(a)<f>j(rW) are the wave functions transformed by 
the symmetry operator R(a) and Oij(a) is the overlap matrix: 

Oyia) = ! tfr^tfircT^faira). (51) 

Finally, for the RPA amplitudes given by non-hermitian matrix one has 

^V)=£^(rV). (52) 

3 

To derive the fields, first we recall that the variation of the non-local 
density, over the wave function reads 

= J]J]<(7'|(7„„|(7)A(r(7)5Vi(rV). (53) 

i era' 

Therefore, variation of the primary density is given by 

Spn'L'v'j'w = y^y^{[K n >L'^v'-, (r >< r ]j'M'^i{rcr)S't/3i(r'cr')} rl=r . (54) 

i aa' 

At this point, operators K n iy mix derivatives acting on the variables r 
and r', cf. Eq. ( 118ft . By using the Wigner-Eckart theorem, we can express 
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them as sums of products of derivatives D mIMl and D' m , rMI , which act on r 
and r', respectively, that is, 

Kn'L'M' L = Kmlm'I' & IMjI' M'^mlK^D ' m , j, M |, (55) 

mlm'I 1 AfjAff 

where the order of derivative is conserved, 

n = m + m', (56) 

and the triangle rule of angular momentum coupling must be obeyed. Nu- 
merical coefficients K^ m , r can be calculated by using methods of symbolic 
programming. At N 3 LO, only 91 coefficients K^w 

are needed, so they can 

easily be precalculated and stored. 

We now can insert expressions (j5"ll) and (155]) into Eq. (jl5l) and remove 
condition r' = r by adding the integral over r' and the function 5(r — r'), 
that is, 

5£ = E [ d 3 r f d 3 r'5(r-r') 

n'L'v'J' ^ ^ 

x ^^[Un>L>viji{r)[K n , v a v ,. a , a ]j,}Q(l) i (ra)5<ilj i (r'(j'). (57) 



This allows us to integrate by parts over r ' and transfer the action of D' m , IIM , 
onto the delta function. With all the angular momenta couplings shown 
explicitly, this gives 

SS = E f d * r f dV E E { D 'l'i'M^(r-r')} 

n'L'v'J' J mlm'I' MjM'j 

x EEEE C^/^-U n , LV] ,,_ M ,(r) 

i era' M' L fi' M' 

x ^mim'i'C IMl \ tM ,pmiMj(yv'a'^'^i{ra)5i)^r'a'). (58) 

The action of D' T mlI i M ^ onto the delta function can be replaced by that of 
(— \) m ' D^, j, M , , and then the integration by parts over r gives, 



58 



E /dv/dV E E <^-0 

n'L'v'J' J J mlm'I' MiM'j 
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X Y Y Y Y C L>M' L v'»> ( V2J'+1 (~ 1 ) m D mTM'U n ' L 'v'J'-M'(r) 
i (tct' M' L ix' M' 

x K^r I ,Cf^, M/i D mIMl a v ^ l . <T/a (l) i (ra)5'4) i (r , a'). (59) 

In the resulting local integral we request that the operator acting on 
4>i(ra) is equal to the mean- field operator, which gives 

Hp) = Y Y Y Y Y c l'm' l v'»> ( vL'+i 

n'L'v'J' mlm'I' Mi Mj M' L /J.' M' 
x ( — 1 ) m K^f m , r C IMi f, M iD m i i> m'j U n 'L' v> J' , - A/' ( V ) <J V > a ' D m j Ml ■ ( 60 ) 

In this form, the mean-field operator is expressed through sums of derivative 
operators standing on both sides of the potentials Un'L'v' j> •■ This form 
can easily be used in the calculation of the matrix elements, because the left 
derivative operator can simply be applied onto the left wave-function through 
the integration by parts. Moreover, potentials U n 'L'v'J',-M' are related to the 
secondary densities by very simple relations f!46p . However, it turns out that 
numerical calculations are much faster if the derivatives appear only on one 
side of the potentials, like it was postulated in Eq. ([7]), that is, 

HP) = Y [ U nWj>(r)[D 

n'L'Gv'] J']o 

n'L'v'J' 

= Y Y ^ V2J'+1 CifM' L v'ii>Un'L>v'J\-M'(r)(Tv'»>Dn>L'M' L - (61) 
n'L'v'J' M' L n' 

3.4- Fields for terms containing additional density dependence 

In the Skyrme functional, a dependence on the isoscalar density is usually 
added to the four terms that are zero order in derivatives, p\ and (r = 0,1), 
that is 

^0,0000 / \ _ ^ya,0000 / \ a 2 / R r>\ 
J oo,ooool r J ~~ ^OO.OOOO \ r ) PoPti l 0Z J 
ma, 0011 / \ _ ^a,0011 / \ ot 2 / ao\ 

1 00,0011 V ) — O 00,0011 \ T ) Po s t- 

We use the extra index a on the coupling constant to distinguish the notation 
from the one used for the density-independent terms. 
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The contribution to the fields is obtained from the variation of the energy, 
which we show here only for the first of the above two terms, namely, 

rpa,0000 / \ o /^-ya,0000 / \ a \ o 2 

oc -Q0,0000 V) _ ° l°00,0000 V) Po) 2 i ^ya.OOOO / \ a °Pt 

5<j>* (t>) d<j>* (r') Pr + G ° ' 0000 [T) Po d<f>* (r') 

= Cd°o o o°o ( r ) [aPo~VX',o + 2p£pA v ] (pi (V) 

which gives the additional contributions 
1 

£ C oo o°o (r) [apT + 2p?pAv] (64) 

T=0 

to t/oooo,o (t')- 

3.5. Rearrangement terms 

For density- dependent terms, the total energy obtained from the proton 
and neutron eigenvalues ef' n and kinetic energies T p,n 

i 

includes the additional rearrangement term 

Err = \{T p + T n )- l -Y, ^ + ^) + / d 3 r H(r). (66) 

For spherical symmetry the terms with density-dependent coupling constants 
in the Skyrme functional only involve the p T densities and for these a straight- 
forward derivation gives 

3.6. Densities 

In order to calculate the potentials in Eq. fH6|) we have to determine all 
secondary densities: 

PmI,nLvJ,J'M'{r) = [An/ PnLvj(r)]j'M' = [An/ [K-nLpv Ol , f* 2 )] j] J'M' ,(68) 
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where D mI is built by coupling gradients Vi + V2 and K nL is built by coupling 
gradients (Vi — V2)/2z. Here and below we understand that r = 7*1 = r 2 is 
set after performing all the differentiations. 

In analogy with Eq. fl55l) . we can split these operators as 

K nL M L = ^2 K?Rr'R> ^2 ^RM^R'M'J^rRM^lm'M'^ (^) 
riir'i?' M R M' R 

D m iMi = (2«) m ^ (— 1) P K™p plpl ^ CpM p p>M' p -DpPMpDpip'M' p '(7Q) 
pPp'P' MpM' p 

where n = r + r' and m = p + p'. Then the density ( I68p reads 

PmI,nLvJ,J'M'(r) = ^ C'lM^Mj C LM'[vu ( 71 ) 

MiMj M L /i 

xfOftm (_r\p' K mI r /M 7 ^(1) D (2) 

x l zi J v ' ^pPp'P' ^PMpP'M'p^pPMp^p'P'M'p 

pPp'P' MpM'p 

K nL V" r iMi ^(1) D (2) , s 

X ^rRr'R' 2-*i RMpR' M' R rRMp r' R' M' R r v ^\ 1' 2/' 

We now introduce coefficients D^ pP , which allow for expressing products 
of derivatives as: 

^rRMp^pPMp — 2-^< 'RM R PM P - u rRpP- u uUM u i 
uUM v 

n( 2 ) n( 2 ) -V s r u ' M u n u ' u ' (7*\ 

^r'R'M'p^p'P'M'p ~ 2.^1 ' R'M'pP'M'p^r'R'p'P'^u'U'M'^ \'°) 



'i'U'M'jj 



where u = r+p and v! = r'+p', that is, u + u' = m+n. These coefficients can 
be calculated by using methods outlined in|A] At N 3 LO, only 91 coefficients 
DrRpp are needed, so they can easily be precalculated and stored. The sum 
of products of four Clebsh-Gordan coefficients can now be recoupled (see 
Eq. 8.7(20) in Ref. Q) as: 

Ey- f-iIMi fiLM L r VM v r U'M' v 

/ j PMpP'M'p RMpR' M' R RMpPMp R' M R P'M' p 
MpM'p M R M' R 

= ^(2/ + 1)(2L + 1)(2U + 1)(2U' + 1) 
( P' P I \ 

X 1 R> R L } ^MuV'M'pLM^IMj ( 74 ) 

wm w U' U W 
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Subsequently, the sum of products of three Clebsh-Gordan coefficients can 
be recoupled (see Eq. 8.7(12) in Ref. [6]) as: 



E r WM w r J'M' piJMj 
° LM L IMi^ IM i JM j Lj LM l v/i 

M j Ml Mi 



■1)^-V(2^ + 1)(2J+1){ L Jt I W j 



• (75) 



The last two remaining Clebsh-Gordan coefficients can be absorbed in the 
following definition of the coupled derivative of the density: 

pl U j^' W {r) = [[D®D%] wPv {r 1 ,r 2 )]j, u . (76) 
_ v r J ' M ' V r WMw n {1) n {2) n (r, r,\ 

~ ^WMwVfi 2-^1 ^UM u U'M[ r - U uUM u - U u'U'MljP'"A' 1> ' 2>l , 

M W fjL M V M' U 

which finally gives 

_ \^ A uUu'U',W n uUu'U'W f„\ ( 77 \ 
PmI,nLvJ,J'M'[r) - ^ A mI,nLvJ,J' PvJ> M> l r J> I' 'J 

uUu'U'W 

where coefficients A^j'^Jjji resu lt from summing up all intrinsic indices: 

auUu'U',W _ (n-\m ST^ \^ ( i W jsml jsnL j^uU Tyu'U' 
A mI,nLvJ,J' ~ 2-^1 2-^i ^ > IX pPp'P'- PL rRr'R'- u rRpP- U r'R'p'P' 

pP p /pi rRr'R' 

( P' P I 

xy/(2I+l)(2L + l)(2U+l)(2U' + 1) I R' R L 

{ U' u w 

x(-l)^V(2^ + l)(2J+l){ j, [ W j). (78) 

At N 3 LO, only 3138 coefficients A^^ l ^ji are needed, so they can easily be 
precalculated and stored. 

4. The N 3 LO potentials, fields, and densities in the spherical har- 
monic oscillator basis 

4-1- Spherical HO basis 

The standard spherical HO wave functions, are given by 

<!>m jm (r,e,<i>,<r) = tf /2 F m (br)e-& )a V Cf\ Y lmi (9,<f>) X i (<r),(79) 



lmi7;m s 
minis z 
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where b is the oscillator constant 



-jp (80) 



and Fj^i(br) = are proportional to the standard Laguerre polynomials [zj]. 

To calculate the secondary densities (176]) . one has to act on the space part 
of the wave function (1791) with the derivative operators D nLML . This leads 



to defining the polynomials m (br) such that 



nLM L Nlmi 1 



^nLA/ L 6 3/2 ^(6r) e -|^V /mi (i 

m+3/2 ^ r ~fcm, /, n (U -i- 



E F ^M A (He^^n mfe (^0), (81) 



A' 



where = + m/. From the Wigner-Eckart theorem, these polynomials 
must have the form: 

^nLM L Nlm l ^ >r ) = ^2k+l ^lmiLM L ^nLNl (^ r ) ' (82) 

and our goal is to determine the set of reduced polynomials F^ LNl (br), in 
terms of which the derivatives of the spherical HO wave functions read 

D nLML <p mjm (rAM = b n+3/2 J2 E VffiCM^ lms 

1 , 



x F* LNl {br)e-2^ Y kmk {9A)x± m {v) 



x {!• i ; )^(He^ (br) v fcmfc (^0) Xlm (a). 

(83) 

Explicit form of F^ LNl {br) can be calculated by using the Wigner-Eckart 
theorem again, namely, Eq. ( 1HT1) must have the form 

D n LM L 4>Nlmi = E 7 2 fe^^LA/ L (^A r 'fcll^ill^) < / , A r 'fcm fc , (84) 
N'km k 
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where 4>Ni mi is the space part of the wave function (1791) . Then we have 
polynomials F* Lm (br) expressed through the Laguerre polynomials as: 

b n F^ LNl (br) = J2^'^ D ^Ni)F N/k (br), (85) 

N' 

where the reduced matrix element can be calculated by considering only one 
matrix element, namely, 

(0Af'fc,m fe =o|-DnL,A/ i =o|0A r «,mi=o) = ^ 2 fc+l ^MLO (4>N'k \ \ D n L 1 1 4>Nl ) ■ (86) 

Note that the Clebsh-Gordan coefficient C^q L0 is not zero for angular mo- 
menta restricted by the parity conservation, (— l) L+z_fc = 1. 

The parity conservation induces specific conditions on the polynomials 
F£ LNl (br). Indeed, by comparing parities of both sides of Eq. flHTT) we see 
that 

F k nLNl {br) = 0, for = -1. (87) 

Equivalently, since for all derivative operators we have (— l) n+L = 1, we see 
that 

F k nLNl {br) = 0, for (-1)*+"* = -1. (88) 

Since polynomials F N[ (br) are real, phases of polynomials F^ LNl {br) are 
fixed by those of the derivative operators fl3T|) and spherical harmonics ((J, 

ijw(M) = (-l)- M Yj,-M(e,<P), (89) 

that is, 

F&mibr) = {-l) l ~ k F* LNl (br) = {-l) n F* LNl {br) = {-l) L F* LNl (br), 

(90) 

where the last two equivalent forms result from Eqs flHTj) and fl88|) . 

The F^Tjft polynomials are calculated using formulas for spherical deriva- 
tives (see p]) combined with recursion relations for derivatives of Laguerre 
polynomials. In this way spherical derivatives of one of the basis functions 
(pNijm ( r , o) = 9ni ( r ) ±Y\ m (9, (p) can be expressed as a sum of functions 

i 

Vi^.-Vi^nj (r)lj ffl = ^ [a{l,n,i 1 ,..,i N ,r)g n _ 1: i(r) 

il,..,ijv=— 1 

+ b (I, n, ii, .., zjy, r) g n \ (r)j Yz_|_i 1+ .._|_i njm+Atl+ .. +AtiV , 
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where the a and b coefficients needed for the different orders were derived 
using symbolic programming. In this way we obtain analytical expressions 
for all derivatives which can then be calculated with good accuracy. 

4-2. Densities in the spherical HO basis 

For any density matrix one can always perform a multipole expansion. 
This strategy fits very well our applications to the spherical HF solutions, 
where only the monopole component of the density matrix is nonzero, and 
to the RPA applications in spherical nuclei, where all multipole excitations 
separate from one another. Therefore, in what follows we consider the density 
matrix of multipolarity J in the form given by the Wigner-Eckart theorem: 



and depending on its reduced matrix elements (4>mj\ \p J \ \<pN'i'f) Then, the 
non-local densities can be expressed in terms of the spherical HO wave func- 
tions ( !79l) as 




(91) 



Pi' M "(n,T 2 )= <t>Nljm{ri,Cr)p[ 



Nljm^N'Vj'm'rN'l'j'm' 



l (r 2 ,a / )(a / \a Vfl \a), 



Nljma 
N'l'j'm'a 1 



(92) 



that is, 




X 




x F N , v {br 2 ) C3 ,' m ^ ,~YvJ&. 



, 2 >0 2 )xi m/ K)- (93) 

2 s 



We can now replace the spin coordinates by the spin projections, 





(94) 
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and we can use the Wigner-Eckart theorem for the Pauli matrices, 

(XW§™.> = (|INI§>, (95) 

with 



(|lko|||> = v^ , (ilkiiii) = -*V6. (96) 

After inserting the nonlocal density (1931) into the expression for local densities 
(I76p . and after acting with derivatives on spherical wave functions, as in 
Eqs. (JED) and ([82}, we obtain 

~uUu'U'W,J"M" , n _ n J'M' ST^ n WM w / ^U'-M'r, 

PvJ'M' {' )— / 4 ^WMwvu / 4 ^UMjjU'M^y L ) 

x y b 3 + u +U ' e - {br ri c l< a\\a v \a) 

* — ' V2 7sm s vfi 

Nljmm s ^ 
N'l'j'm'm's 

x E 



x 72^^wj»A/»(^illP J ''ll^'^') ( 97 ) 



7 O ' 



where, by using the phase convention D j v , M i = (— 1) D K J V , _ M , , we 
have introduced the complex conjugation into the derivative operator (!3~T1) . 

After a lengthy but straightforward derivation presented in IEJ we obtain 
the following result: 

~uUu'U'W,J"M" ( \ ~uUu'U"W,J" \ -(br) 2 ^TM T v M j,\ /r^ 

P„J'M' ( r ) = 2^ P vTJ' ( br ) e C JiM>J"M» Y TM T {0,9), (98) 



T 



for the radial form factors p^j /U W ' J (br) given by: 



~uUu'U'W,J" (u r \_( i N 1+t. / 1 1 1 1 1 1 \ / (2 J'+l) \p (_,)j'+J F >ljkU,WJ" 

Nljk 
N'l'j'k' 

xfe 3+u+u >^(^)(^||/5 J 1|0iv^)^W^(^) ) (99) 
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where 



R ljkU,WJ" _ ( _-\ \k( _-\ \T ( i\U-W r T0 / (2fc+l)(2fc'+l)(2iy+l)(2j+l)(2j'+l) 
n l'j'k'U',vTJ' — \ L ) \ L ) \ l ) ^fcOfc'Oy (2T+1) 



£(-lf'(2T' + l) 



x 

T' 

x { I \ \ U V k> U> \ { V T T j, (100) 
I V T' ) [ T' T W 

In view of the fact that only the coefficients B^^l^/^j, for (— iy+ k + u = l 
and (— iy'+ k '+ u ' = l are required in Eq. we may replace them by 

coefficients Bj^'JfJj,: 



BljkU,WJ" _/_ 1 Vr ,T0 / (2fc+l)(2fc'+l)(2W+l)(2j+l)(2j'+l) 
n l'j>k'U',vTJ> — y l ) ^fcOfc'Oy (2T+1) 



$>T' + 1) 



x 

T' 

/' z r \ \ w u w " " 



where we have used symmetry properties of 9j symbols under the trans- 
position of rows and columns and transposition with respect to the main 
diagonal. 

Finally, all secondary densities of Eq. (ITT)) can now be calculated in terms 
of one compact expression: 

f&LjsM'tr) = £ Pl J i M A hr y~ {br?c VM T 'j^ 

where the radial form factors read 

Pm[,nLvj,j'(br) = £ A^Jl^j, p™Tj> U W ' J (br). (103) 

uUu'U'W 

4-3. Potentials in the spherical HO basis 

In the spherical basis, the secondary densities to be used in Eq. (|46p have 
the form given in Eq. (I102p . Therefore, the potentials in Eq. (1601) acquire the 
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form 

^n'LV , j'M'{ r ) = ^Iw;'(^) e ^ Cj>ai>j>'m''Ytm t (8,4>), 

(104) 

where the radial form factors read: 

^WJ'(^) = \pmI,nLvJ + ( — 1) J ^ml^n'L'v'J'j PmI,nLvJ,J r (pT'\105 S ) 

mI,nLvJ 

Similarly, potentials U£' L ¥ V " J/M ,(r) m Eq. (TBTT) read 

u£'£!j.M'tr) = £ ^i;^'(^)e-^ 2 c™v M „FTM T (^ 0), (106) 

TAf T 

and the radial from factors U^, J L \ vlJ ,{br) can be calculated in the complete 
analogy to the results outlined in Section [21 see Eq. (jSj), namely, 

Un'L'v'J'(br) = CP a Xa] a -,n'L'v' J J'Pm J I,nLvJ,J'(br). (107) 

aafi\ml ,nLv J 

4-4- Matrix elements of the Hamiltonian liUty) in the spherical HO basis 

In the spherical HO basis of Eq. (1791) . we can calculate the matrix elements 
of the single-particle Hamiltonian (|60|) in the following way: 

{N'l' 3 'm]\KP J '' M '')\Nl 3 m J )= £ £ ^ ^ U f dfi^ 

n'L'v'J' mlm'V MrM'j M' L fi' a 1 a 

E r J'M' {-\Y'- M ' ( i\m' Tsn'L' r<L'M' L , n , , 
u L'M[«y V2J'+1 ^ ' JX mIrn'I' Ly IM I I'M' I \ <J \ a v' u' \U / 

M' 

4>*N>l>j>m' 3 ( r , 6 > 0> °"0 Dm'I'M'finWJ'-M'i^DmlM! 4>Nlj mj (r, 0, 0, <j), (108) 

where h(p J " M ") denotes the field calculated for the one-multipole density 
matrix p J M (19~T1) . The integration by parts now gives 

n'L'v'J' mlm'V M 7 M| M{y ct'ct 

E^J'M' (_i)J'-m' ,„,/ „/!) ^i'-Mz, / / 1 I \ / iW+Z'-M' 
°L'JI/[d>' v/2J'+1 > 1X mIm'I' U IM I I'M' I \ <J \ a V fi' \V / K~ 1 ) 

M' 

Un"vv'J'-M'( r ) {Dm'I'-M'^N'l'j'm'^r, 8, (f), (D mIMl (f) Nljmj (r, 6, 0, cr)) , 

(109) 
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where we have used the hermitian-conjugation property (152"j) of the differen- 
tial operators. 

By inserting potentials (11041) and derivatives of spherical wavefunctions 
f l83|) we have: 

(Nrj'm'/h(p J, ' M '')\Nljrn 3 }= £ £ £ £ / r 2 d r /" dfi £ 

n'L'v'J' mlmT MiM'j M' L /jf <y' a 

ECiJ'M' (-i)J'-M' , nra '„„' L ' ^yi'M^ , I w U rn'+I'-M' 

^L'M' L v'ii' ^ 2 J'+1 ^~ > IS -mIm'I'^IM I PM' J \ a \ a v' n' \ a / K~ 1 ) 

W 
TAf T 

v m'+3/2 \^ 1 n k ' m 'k r j ' m 'i 

XD 1^1^ VW+l U Vrn[I>,-M^ Vm ,l 

k'm' k m',m' s 

x^™(&Oe^ (M V fc * m ,(^0)xi m/ (a') 

1 



x^ m (6r)e^ W WM)xl m (a), (HO) 

The angular part can be integrated explicitly, by using the multiplication law 
of spherical harmonics (Eq. 5.6(9) in Ref. |6j), and summations over a and 
o' can be performed as in Eqs. (|94l)-(!96l). This gives 

{N'l'j'm' ] \h{p J '' M '')\Nljm J )= £ £ ^ ^ f r 2 dr 

n'L'v'J 1 mlm'I' MtM', M' t a' 



E n J'M' (-!)■'' ~ M ' i .sm'r^n'L' n L' M' L , - 
u L'M' L v'n' s /2Ji+\ \ L ) Jx mIm'r'^IM I I'M' I \ . 



m'+I'-M'j 



EV V" V" i r 2 m s /i|i_ II i\ / (2r+i)(2fc+i) r yfc , o 
2^ 2^ V5°I_ _/..,\2ll <7 « / ll2/V 47r(2fc'+l) TOM TMrfa|j 



AT 

I™' 

^ \ * \ ■ \ * \ ■ \ * i 

Z — <• Z — / Z — i L — / Z — / J- m 

TMt k'm k mjm' s fem fc mim s ■« 

x & 3+n >^^,(^)^iv J ,(6r)F^ 7V/ (6r)e- 2 (^) 2 

^J',-M'J"M" VW +I U l'm' l I',-M^ l , m , l m , ^S+I U »m,/M^ Imj i m# . U^J 

where we have used condition (f5B"j) . 
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We may now proceed by calculating the reduced matrix element of the 
field: 



E Tmrf^"^ {N ' l ' fm 'i Wf'nmrui). (112) 



rrij m'j Mj 



>N>i 

X 



Again, after a lengthy but straightforward derivation presented in we 
obtain the following result: 

n'A\h I "{~P J " M ")\\<t>Ni j )= h»A-l) J " 

E E ^/'(-ir'(lii^iil)v^E(- 1 ) r ( 2T+1 ) 

n'L'v'J' mlm'P ' Tkk' 

x J rMre- 2 ^ 2 b* +n 'F* P Mbr)U^ 

R ljkI,L'J" , m \ 

where we see the same numerical coefficients ( 11001) that already appeared in 
Eq. pi. 



^.5. Matrix elements of the Hamiltonian |7|) in the spherical HO basis 

For the Hamiltonian in the form given by Eq. ([7]), a similar derivation 
gives the matrix elements in the spherical HO basis that can be written as: 



i. j ,\\h I "{p J " M ")\\<i>N lj )=5i»j» 



x / r 2 



x 

v'J'T 



drt 3+ "'e- 2W2 C(6r) E F ^m(br) 

n'L'k 

Y,VS2*A*)(l\Mi)h$$Z J '' (H4) 



with the reduced form factors U^', v ,j,(br) given in Eq. (11071) . and 

,ljkL',J'J" _ 1 ( i \J"+l+v'+l'+L' 

v /(2k + 1) (2f + 1) (2j + I) 

' l i'j',v'T - i; v^FTT) 

x (-l) T (2T + l)Cl' ° T0 

f 3 3 «' 

x V(-i) T (2T' + iW y j j" 

T' [ V I V 
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r r i V \ j r t f \ . 
x \ k t v } \ v v> r j • 1 b) 

4-6. The total potential energy in the spherical symmetry 

The energy related to one term in the energy density (1381) reads 

cn'L'v'J' _ f A 3 rr n'L'Jj> ST j-l) J '- AI ' ( A™2 \( n J " M " \ 

^mI,nLvJ I ^mI,nLvJ ^2J'+1 \"00,n'L'v'J',J'M'J y^ml ,nLvJ,J' ,-M' J ' 

(H6) 

We carry out the derivation for arbitrary values of multipole components, 
that is for J'^M'^ ^ J'{M'{ ^ 00, while the standard HF total energy in spher- 
ical symmetry corresponds to the particular case of J'^M'^ = J"M" = 00. In 
the general case, the result corresponds to the energy matrix element between 
two multipole components, and allows one to calculated the total energy for 
the density matrix being a linear combination of multipole components, that 
is, for a deformed state. 

After using the definition of reduced density in Eq. (I102p . the energy 
(11161) becomes equal to: 

cn'L'v'J' _ f A 3 m n n'L'v'J< „-2(br) 2 \^ (-1) J '- M ' 
^mI,nLvJ I ml ,nLvJ^ / j ^/2J'+1 

J M> 

X P^%L'v'J',J'( hr ) CT J'M^J''M^T'M' T {0,4>) 



T'M' T 



K ^ PmI,nLvJ,J> (br)C™J M , j,, M „Y TM/j \n.< 



TM T 



By using the multiplication theorem for the spherical harmonics, one obtains 
the expression 

cn'L'v'J' _ f ,3 rC n'L'v'J' -2(brf ST h T ' J 2 (hr\ ?? J " (hr) 

C mI,nLvJ — / U T °m;,nL»; e / , rQQ.n'L'v'J'.J' \ ' ) Pml.nLvJJ' \ ' ) 

J t"-tm 



TT'Q 



X 



(2T + 1) (2T> + 1) ^ (-!)■'' -m' T 'M' T 

47r ( 2 Q +1) 2^ V2T+1 W'M',J»M» 

v ^ ' M'M' T M T M Q 

x Cj, J M ,jn M „C TM ® TIM ^YQ MQ (6,<p) , (H7) 
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which, for density-independent coupling constants, can be integrated over 6 
and (j). After summing up the Clebsh-Gordan coefficients, one obtains: 

cn'L'v'J' _ j-l) J l- M " x r f 2 j „ r n'L'v' J' -2(br) 2 

t m I,nLvJ ~ V2 J'+l(2J{'+l) d J2Ji d M^-K' J r dr ^mI,nLvJ e 

x X)(- 1 ) T (^ +1 )^w^(^)&t,J^ (^)- ( 118 ) 

T 

One thus obtains the correct result that for density-independent coupling 
constants, the potential energy is diagonal in different multipoles. 

When the coupling constants do depend on density, the angular integral 
cannot be performed, and a similar derivation gives the general result: 

cn'L'v'J' _ f d 3 nn'L'v'J' -2(br) 2 f J'i /.n ~TJ[' /,x 

C mI,nLvJ — / a ' ^ml^LvjZ / , PoO,n' L'v' J' ,J' \ ' ) Pral \nLvJJ' \ ' ) 

J TT'Q 

Q0 (2T + 1) (gr + 1) 



'TO,™ /47i 



X 



( J' J' 

E°WW T T Q \Y QMQ (9A). (119) 



For the spherical case, one has J" = J'l = 0, which implies T = J 1 , and 
both expressions f 1 1 1 8 j) and (If 19[) reduce to: 

cn'L'v'J' _ /_-. \ J' /q 7/ _i_ 1 fjj. nn'L'v'J' p -2(br) 2 
C mI,nLvJ — \ l ) V ZJ + 1 / r Qr °raf,iiW e 

X P00°n'L'v'J',J' (^ r ) Prrd,nLvJ,J' (^ r ) ■ (120) 

For density- dependent terms, such as those in Eqs. (162|) and (|63|) . the 
total energy reads 

cce.n'L'v'J' I i\J' / 7/ ; T^ya,n'L'v'J' I 2j -2(M 2 ol ( ~l \ 

£ m i,nLvj = (-l) J V2J' + lC m ; inW Jrdre { > p Q (br) 

X P00°n'L'v'J',J' (^ r ) Pird,nLvJ,J' (^ r ) 5 (121) 

where the density po (& r ) depends on the reduced density as: 



Po(H =^00,0000,00(^)^7^- ( 122 ) 
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4-7. Direct Coulomb energy in the spherical symmetry 

The integral for the direct Coulomb energy in spherical symmetry has 
a integrand which has a discontinuous first derivative. In order to get an 
integral which is easier to calculate using Gauss-Hermite integration we use 
the 'Vautherin trick' [8j] 



I dVdVp(r) , 1 , p(r') 
J \r — r \ 



2 „ 
(4tt)V 
6 

47re 2 



dr dr' 



(r + r 



2\/3 
x dr 



dr dr'poo,oooo,o (r) re 



-{brf 



rr'p (r) A'p (r') 



fr + r 



_/|3 



r'p 2 o,oooo,o (r') e^') 2 



Aire 2 
2^3 



Poo.oooo.o (7) 



(123) 



which gives a smoother integrand. However in order to perform the linear- 
response calculations, where expressions for non-spherical multipoles are 
needed, we also developed a different method which will be presented else- 
where [9|. 



4-8. Exchange Coulomb energy 

For the exchange Coulomb energy we use the Slater approximation 



E P 



3 /3 



4 \7T 

3 /3 



1/3 



4 V^r 

3 4 / 3 e 2 



1/3 



e 2 47r / r 2 drp p 



.4/3 



e 3 



>r) 2 



(4) 



2/3 



6 (br) 4/3 ei (fer)2 

^P00,0000,0 v U ' ) e 



(4tt) 



The matrix elements are the same as with density-dependent forces (using 

a = -2/3 and C^SSSo = ~f (f) V3e2 ) but should onl y be added to the 
proton matrix elements. 
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4-9. Numerical integration 

Numerical Gauss-Hermite integration is used to calculate the radial inte- 
grals occurring in the expressions for matrix elements [Eqs. (11131) and (I114p 
and total energies [Eq. (11201) ]. This kind of integration is for integrals of 
the form e~ x f (x) dx, and in order to obtain this form our integrals are 
transformed by using x' = (V2b) x. The integrals can then be written as: 

N/2 
" i=l 

Where f sym (x) = f (x) 9 (x) + f (—x) 6 (— x) was used in an intermediate 
step. To reduce the number of grid points by half to iVg r id = N/2 it was 
also used that for Gauss-Hermite integration, the weight functions Wi and 
grid points x" = = bri are symmetric about the origin. The integrals 
for matrix elements and total energies of most terms become exact when 
A/grid = iVo + 2, where Nq denotes the maximum HO shell included in the 
basis. But in general more points are needed when the integrand cannot be 
expressed as a product of four basis states, e.g., in the case for the Coulomb 
interaction and also for the density- dependent terms. 

5. Overview of the software structure 

In order to get an idea of the software structure we list some of the main 
calls performed by the program. In the code hosphe (vl.00), there are also 
deeper level calls that are not listed here. Depending on the input data, some 
of the calls listed below may or may not be performed, but all are anyway 
included in the list. We specify the subroutine calls by first writing the name 
of the module followed by the name of the subroutine. 

• hfmain 

Driver program which reads input and starts the run by calling the hf 
routine. 
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- hf 

Main program which contains the iteration loop. Starts by calling 
routines which initializes different tables. 

* RecoupDK: define _cfk 

* RecoupDK: define_cfd 

* RecoupA: define_cfa 

* RecoupA: define_cfa2 

* Coupling: restrict_densities 

* Coupling: get_cc_restrict 

* Coupling: makeJtab 

* Skyrme: get_skyrme_cc 

* Coupling: read_carlsson_from_file 

* Coupling: cc_perl_to_n31o 

* Coupling: set_dependent_cc 

* Coupling: cc_n31o_to_perl 

* grid: set_grid 

* qnotra: set_sps 

* qnotra: qpbase 

* mflib: Read_hf_initial_density_matrix 

* mflib: set_onedet 

* Ho_Derivatives: store_derivatives_of_basis 

* Ho_Derivatives: read_derivatives_of_basis 

* TypeDefinitions: inverse_derivative_densities 

* F2U: Define_F2U 

* mflib: read_density_matrix 

Next follows the main iteration loop. 

* Derivative_densities: sum_derivative_densities 

* hmatrix_Ud: hamiltonian_matrix 

* mflib: diagonalizeJif 

* mflib: hf2ho 

* density _energy: energy _from_densities 

* density _energy: energy _from_densities_densdep 
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* density _energy: energy _from_coulomb 

* mflib: hLenergy 

Finally there is some post-processing of the results. 

* density _energy: plot_energy_terms2 

* sorting: sortrx 

* mflib: store_density_matrix 

6. Description of the individual software components 

The code hosphe (vl.OO) is separated into different files with each file 
containing a module. Modules may contain module parameters and collec- 
tions of subroutines and functions. Below follows a list of the files defining 
the code where the *'s should be replaced by the latest version number. 

The modules which are used mainly in the initialization phase are: 

• hfmain_dist.v*.f90 

Driver routine which reads the input and calls the hf subroutine. 

• hfmain_Parameters.v*.f90 

Allocates parameters used as input to hf. 

• qnotra.v*.o 

Defines quantum numbers for basis states and quantum numbers for 
the reduced density matrix. 

• skyrme.v*.f90 

Contains coupling-constant sets for various standard Skyrme EDF's 

• coupling. v*.f90 

Assigns quantum numbers to the coupling constants. 

• TypeDefinitions.v*.f90 

Contains list of the allowed quantum numbers for the derivative oper- 
ators, local densities, and derivative densities. 

• RecoupKlist.v*.f90 

Contains one subroutine which reorders the list of K coefficients. 
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• RecoupDK.v*.f90 

Contains values and quantum numbers of the D and K coefficients. 

• RecoupA.v*.f90 

Contains values and quantum numbers of the A coefficients which are 
used to recouple derivative densities. 

• F2U.v*.f90 

Contains the F2U coefficients. These are used to obtain fields for each 
term in the functional and are almost the same (see Table [7]) as the \ 
coefficients in Eq. (|S]). 

• Ho_Derivatives.v*.f90 

Contains routines to calculate derivatives of the HO basis functions. 

• grid.v*.f90 

Generates Gauss-Hermite points and weights. 

The modules which are used mainly in the iteration loop are: 

• hf.v*.f90 

The main subroutine containing the iteration loop. 

• Derivative_densities.v*.f90 

Generates the reduced derivative densities as functions of the grid 
points given a reduced density matrix as input. 

• hmatrix_Ud.v*.f90 

Determines the matrix elements by using for the fields the form in 
Eq. (J7|). Calculates the reduced matrix elements given reduced den- 
sities and reduced matrix element quantum numbers as input. Also 
calculates the matrix elements of the density- dependent terms and of 
the Coulomb interaction. 

• matrix.v*.f90 

Determines the matrix elements by using for the fields the form in 
Eq. (J60~]) . Calculates reduced matrix elements given reduced densities 
and reduced matrix element quantum numbers as input. Used mainly 
to test the accuracy of the matrix elements from hmatrix_Ud.v*.f90. 
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• density _energy.v*.f90 

Calculates integrals by using the Gauss-Hermite integration. Contains 
subroutines to calculate integrals over the terms in the functional to 
obtain the energy. Also calculates Coulomb energy integrals and radii. 

• coulomb. v*.f90 

Calculates Coulomb direct and exchange matrix elements. Used mainly 
to test the accuracy of the matrix elements from hmatrix_Ud.v*.f90. 

• densitydep.v*.f90 

Contains density-dependent factors that multiply each density-depen- 
dent coupling constants in the mean field potentials. 

There are also modules containing helper routines: 

• sorting. v*.f90 

Contains routines used to sort data. 

• geometric. v*.f90 

Contains a collection of routines for the Clebsh-Gordon coefficients, 3, 
6, and 9j symbols, Gamma functions etc. 

• mflib.v*.f90 

Contains a collection of auxiliary routines such as the diagonalization 
routines and various transformations. 

7. Differences between the notation in the code hosphe (vl.OO) and 
in the text 

In some cases, the notation used in this study differs from that used in 
the code hosphe (vl.OO). For example, to store real objects in the code, we 
have removed complex phases from some coefficients. The transformation 
rules for some useful objects are defined in Table H 
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Notation in text 


Definition in the code 


Stored in module 


P;dS 


(i) n '> (-i) n ^ F2U%Coeff 


F2U_table 


r nLNl 


H) n F_val 


Ho_Derivatives 


D n'L' 


D_REC%Coeff 


RecoupDK 


K n'L> 


(2iy n ' (-l)- m 'K_REC%Cocff 


RecoupDK 


auUu'U',W 
ml ,nLvJ,J' 


(2i)-™A_REC%Coeff 


RecoupA 


u 7 


\ (if ILfkn 


hmatrix_Ud 



Table 1: Table relating notation introduced in the text with corresponding objects defined 
in the code hosphe (vl.OO). 

8. Description of input data 

Input is given to the code hosphe (vl.OO) by using the FORTRAN 
"namelist" statement. In this way, the variables specified in the input have 
their values assigned, while those not specified in the input retain their pre- 
defined default values. The variables that can be specified in the input are 
listed below. 

1. noscmax 

Maximum main HO quantum number N used in the HO basis. The 
code hosphe (vl.OO) currently supports values of N up to 70. 

2. order max 

Maximum order in derivatives used. The possible choices are 0,2,4 and 
6. 

3. ngrid 

Number of Gauss-Hermite grid points. The code HOSPHE (vl.OO) cur- 
rently supports up to 85 grid points. If a negative value is given, the 
code uses N grid = N + 2 + 10 points, that is, 10 more points than is 
needed for the most terms to become exact. However, for the Coulomb 
and density- dependent terms to converge with high precision, one may 
need more grid points. 

4. intera 

Name of the Skyrme functional to be used. At present, supported ver- 
sions are: "SLY4" , "SLY5" , "SKM*", "SKP", "SIII", and "FILE". If 
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the name "FILE" is specified, the coupling constants in the spherical 
notation are read from file "cc.inp". An example file of this type is 
included in the distribution. 

5. AZ, AN 

Number of protons and neutrons. 

6. hbarom 

Oscillator frequency fno in MeV. If a negative value is given, it is cal- 
culated as hw = 1.2 x 41A -1 / 3 , where A is the number of nucleons (see 
above) . 

7. boscil 

Oscillator constant b in fm _1 . If a negative value is given, it is cal- 
culated from Eq. (180]) by using the value of Jtlj (see above) and the 
nucleon mass m being the average of the neutron and proton masses. 

8. icm 

Center of mass correction. For icm=0, no correction is used, and for 
icm=l, the code hosphe (vl.OO) uses the one-body center of mass cor- 
rection (E c . m . ~-i(T)). 

9. icoudir 

For icoudir=0, no direct Coulomb term is included, and for icoudir=— 1, 
the code hosphe (vl.OO) calculates the direct Coulomb energy by us- 
ing the Vautherin method, see Section 14.71 

10. icouex 

For icouex=0, the Coulomb exchange term is not included, and for 
icouex=— 1, the code hosphe (vl.OO) calculates the Coulomb exchange 
energy by using the Slater approximation, see Section 14.81 

11. itermax 

Maximum number of iterations allowed before aborting. 

12. epsilon 

Accuracy parameter. The iterations are stopped when the ground-state 
energies calculated by using the EDF and HF expressions differ by less 
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than epsilon and every HF single-particle energy changes less than ep- 
silon between two iterations. 

13. alpha 

Mixing parameter a to slow-down/accelerate the iteration convergence. 
It mixes the density matrix from the current p\ and previous p itera- 
tions, so that the new density matrix is obtained asp = a*p 1 + (i — 
a) * p 

14. keta_J 

Turns the Skyrme tensor coupling constants ON/OFF (1/0). 

15. restart 

If 1 or 2, the code hosphe (vl.00) attempts to read the density matrix 
from the file named as in the following example: 

densities_050_082.rec for (AZ,AN) = (50,82) and the density matrix is 
also automatically stored to the same file when the iterations are fin- 
ished. For restart = no restart of iterations are attempted. For 
restart=2 the stored density matrix is used even though it may come 
from a calculation with a different number of oscillator shells. For 
restart=l it is used only if the number of oscillator shells is the same. 

16. Flag_read_ini_dm 

If .true., the code hosphe (vl.00) attempts to read the initially oc- 
cupied levels from the file "occ_orbs.inp". The first line should say 
N cng , N m where all main oscillator shells up to N m are automatically 
filled and iV cng denotes the number of changes in occupation with re- 
spect to this initial filling. Then follows one line per change, each 
change being specified as N,l, 2j, I occ where N, I, 2j denotes quantum 
numbers of HO j-shells and J occ specifies weather the shell should be 
occupied (I occ = 1) or not (7 occ = 0). This is repeated twice to define 
shell occupancies for protons and then for neutrons. An example file 
of this type is included in the distribution. 

17. verbose 

Verbose is an integer which specifies the amount of output produced 
by the code hosphe (vl.00) during a run. Verbose = is the standard 
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which gives a minimum of output and higher values leads to more 
information being printed to the screen. 



9. Installation instructions 

Most of the source code of hosphe (vl.OO) is written in Fortran 95. 
The code consists of several modules that are linked together by using the 
standard Makefile (see [HI). In general, no special compiler options are 
needed, except for the optimization flags or check flags for tests. These 
options are set in the accompanying Makefile, which is now set to work with 
the G95 compiler that is a free and open-source Fortran 95 compiler (see 

H). 

The code hosphe (vl.OO) also make use of the LAPACK and BLAS 



linear-algebra libraries that can be obtained from [12l . Il3 . 



10. Test run description 

The simplest way to run the code HOSPHE (vl.OO) consists in providing 
the namelist input data in the following form: 

./hosphe « end > pb208.n50.out 
feinput AN=126,AZ=82,noscmax=50, 

icm=l , icoudir=-l , icouex=-l ,keta_J=l , intera="SLy4" ,ordermax=2 , 

epsilon=l . e-7 , itermax=1000 , 

Flag_read_ini_dm = .false. , restart = 0, 

alpha=0 . 65,ngrid=-80,boscil=-2,hbarom=-l/ 

end 

The script above is provided in the distribution file of the code hosphe 
(vl.OO). By executing the script, one obtains the output file "pb208.n50.out", 
which is also provided in the distribution file. The main section of this file, 
which gives the total energies in 208 Pb calculated for the maximum HO shell 



included in the basis of N = 50 and SLy4 Skyrme functional [14|, reads 



Kin.prot Kin. neut . Tot. kin. 

1337.059947 2529.116266 3866.176214 

T=0 Skyrme T=l Skyrme Tot. Skyrme 

-6405.081099 106.598348 -6298.482751 
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Energy 
-1635.692396 



HF Energy 
1635.692396 



Rearr . ene 
1221.821085 



Cou. tot. 



Cou. dir. 



Cou. exc. 



796.614142 827.882912 -31.268770 

In Figs. [1] and El we present results of similar calculations performed for 
208 Pb and the HO bases of Nq = 10-70. Fig. [1] shows the convergence of 
the total energy in function of Nq. It turns out that the energy converges 
exponentially to the limiting value of Eq, namely, 



However, as shown in the two panels of Fig. [IJ two different values of Eq and 
a are obtained for the regions of N Q below and above iVo = 38. A rather rapid 
convergence (a = 0.172(3)) to Eq = —1635.719, which is seen below iVo = 38, 
is followed by a slower convergence (a = 0.1068(8)) to Eq = —1635.69405. 
Since the least-square fit of the limiting values Eq is ill-conditioned, no error 
estimates can be obtained for them. 

By considering convergence patterns in a few more cases for different 
options and nuclei, we found that the trend with two different slopes is not a 
general feature. In the few cases we looked at, we found that it is only above 
40-50 shells that the rate seems to stabilize to an exponential convergence. 
These results show that, in general, its not possible to find the extrapolated 
limit of energy just by calculating only a few points of the curve for some 
small numbers of shells. 

Fig. El shows the dependence of CPU times on N , obtained on the AMD 
Opteron Processor 2352 running at 2100 MHz clock speed. First, one can 
see that the spherical-basis code hosphe (vl.00) is, of course, orders of 
magnitude faster that the 3D code hfodd (v2.40h) [ll]. For N = 36, 
the former needs only 20 sec of CPU time while the latter needs 250 000 sec. 
Second, for both codes, the dependencies on N are clearly given by power 
lows indicated in the figure. Strangely enough, these power lows are different 
for calculations performed below and above iVo = 20. At the moment, no 
explanation for such a timing pattern could be found. 



E(N ) = E + E 1 exp(-aiVo). 



(124) 
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Figure 1: Total energy of 208 Pb as a function of the maximum HO shell included in the 
basis No. Upper and lower panels show in the logarithmic scale differences E(Nq) — Eg for 
two different values of Eq. Large full dots and small empty circles give results calculated 
by using the codes hosphe (vl.OO) and HFODD (v2.40h) [15|], respectively. Up to N — 36, 
where the HFODD calculations could have been performed, one sees a perfect agreement 
between the results given by the two codes. Solid lines give results of the exponential-decay 
fits. 
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Figure 2: The hosphe (vl.OO) CPU times required for calculations performed for 208 Pb 
and the standard Skyrme functional SLy4, shown as functions of the maximum HO shell 
included in the basis No. Circles and squares show results obtained by using the codes 
HOSPHE (vl.OO) and hfodd (v2.40h) respectively. The doubly logarithmic scale in 
the Figure, shows that these times scale as different powers of No, as indicated in the 
figure. 
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11. Summary 



We have developed a computer code that is able to solve the self-consistent 
equations resulting from the use of the generalized N 3 LO energy density 
functional introduced in Ref. For different terms of the functional, the 
resulting mean-field expressions were systematically formulated by using the 
spherical-tensor representation, which was also introduced in Ref. Q. In 
Section [3] we presented the general equations, needed for any choice of basis, 
and in SectionH]we developed expressions needed to implement the formalism 
in the spherical harmonic oscillator basis. 

The present version of the program is based on the spherical harmonic 
oscillator basis and can be used to calculate, for example, ground state prop- 
erties of spherical nuclei using standard Skyrme forces. For this purpose, it 
has some advantage compared to other codes, because of its fast execution 
time and stability. 

The motivation to construct the code, and its main new feature, is the 
possibility to test the influence of higher order derivative terms included in 
the mean field. The values of the higher order coupling constants can be cal- 
culated from the density-matrix expansions starting from other interactions 
or by direct fits to experimental data, and both approaches are presently 
being explored. The formalism was derived in a general way, so as to include 
non-spherical multipoles, and thus it can also be used for linear-response cal- 
culations 16( , which go beyond the mean field. Extensions to include pairing, 
better treatment of the continuum, and statically deformed densities are also 
being developed. 
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A. Calculation of coefficients D™% pP (T72j) 

Coefficients D^ pP ( IT2"j) can easily be calculated by using the fact that 
they exactly correspond to the multiplication of tensors X built from the 
position vector x. 

X rR M R X p p M p = C^M U R PMp D rRpP X uUMu- (125) 

uUMu 
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Tensors X are simply related to spherical harmonics (see Eq. 5.2(40) in 
Ref. |): 

p-p 



X p pm p = 2 r P y/j0^Y PMp (9,<j>). (126) 



Therefore, multiplication law of spherical harmonics (Eq. 5.6(9) in Ref. |6j) 
gives 

R-r 

rUMR- v l>l'Mr y~ y °J 2 r ' V (2A+1)!!' 

P-P 

/ (2R+1)(2P+1) 
V 4tt(2[/+1) 

UM V 
u-U 



X r RM a X v pM B — ( _ 



v r UM v 

x ^ROPO^RMrPMp 



and thus the required coefficients read 

R+P-U 

n uU _ r /^/q^ 2 / i?,!P! (21/-1)!! ^t/Q /inM 

^rRpP — °u,r+p I VOI y (7! (2P-1)!!(2P-1)!! POPO' l ±zc 7 

and are independent of wp. 

B. Derivation of Eq. ( BSD 

We begin by a summation of the product of four Clebsh-Gordan coeffi- 
cients, that appear in Eq. ( 1971) . that is, 

i , 



U l °, 1 ( ~ y j'm'J"M" ( ~\ / ,1 , 

mm s rn'm' g 2 2 2 

= E (-i)^>/^c fB r M i ,(-1)^^^ 

— , V T 7Tm s ^,-m^ V ^ 

mm s m'm' s z z z 

x (-1)^' >JMPf^!U-^ ^ c tU, m > 



mm s m'm' 



mi 
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V C T ' M ' T C T ' M ' T \ x - x - v \ (129) 

/ , Z'mJZ,— m; v,—/j,J",—M" \ 2 2 ( ' 



T'M' T 



see Eq. 8.7(20) in Ref. [6(. 

The product of spherical harmonics reads 



x CCu Y ™ T (^ ( 13 °) 



see Eqs. 5.4(1) and 5.6(9) in Ref. [6j. This gives us another sum of products 
of four Clebsh-Gordan coefficients to sum up: 



%,-m\(~tkm k pk'm' h „T' ' M' T ^TM T 

ImiUMu l'm[U' l'm[l,—mi km^k' ,—m' k 



mim\m k m! k 



E( — '\\ m 'k- rn 'l( — ~\\ l ~ m l I 2k +l (~iUMu 
\ L > \ L ) y 2U+X^km k l-mi 

m,im' l m k m' k 



vf-lV'- m l r^+Lr U ' M 'v r T ' M ' T r TM T 

V / y 2U'+1 I'mlk'-m'^l'm'J-mi km k k'-m' k 

(_l)^-^(_l)I(_l)I'yEJ^(_l)I'^(_l)^-T 



(y u ' M u n T ' M T rjTM-r 
/ j km k l,—mi k' ,—m'ul'm', l'm',l,—mi k' ,—m' k km k 

mi m 'i m k m ' k 

I i\ML-M'/ i\l+k-U'-T 1 2k+l / 2fc'+l 

V" W/Mp r U'M' v r TM T r T'M' T 
X ^km k lmi^k'm' k l'm[ U k'm k km k I'm'jlrni 

mim'im k m'u 

(-l) M T- M v(-l) l+k - u '- T y/(2k + l)(2k' + 1)(2T + 1)(2T + 1) 



w'ml r t w 



Mr 



see Eq. 8.7(20) in Ref. [6j. 

We can now perform the summation over M' v and Mjj, which gives the 
factor (— l) u+u ~ w 5ww'^m w m' and allows for a summation over W and 
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M' w . After inserting all these results into Eq. (157)1 . we obtain 

~uUu'U'W,J"M" ( x ^J'M' flW' ft3+u+u' -(br) 2 1 /111 Mix 

X E 72fa F «W/(^)72fTT(^'IIP J ''ll^'') 

k' 

x (-l) J ' +J - ^+1 -^v / 2JTT^/2JTTv / 2J 7 TT 

f J / J" 

A °«,-fiJ",-W" S 2 2 y 

t'm' t \ i v r 



E / (2fc+l)(2fc'+l) ^TO y- //} 
Y 47r(2r+l) X ^O/c'O^Mt^, 



X 



-l) M ^ +fc - £/ '- T 7(2A; + l){2k' + 1)(2T' + 1)(2T + 1) 
( I k U } 
xC^%m>{ ^ * U> \(-l) u+u '- w . (132) 

( r t w J 

Now we have to sum up products of three Clebsh-Gordan coefficients: 

El _T\M' r <J'M' n T'M' T r WM w 
\ L ) Ly WM w v l i Ly v~iiJ"-M"'^TM T T>M! r 

M w M' T ii 



E, 1 ■s.M'^J'M' / \u+A« / 2T'+1 nJ"M" 

v - 1 / Y 2T +! WM W T',-M' T \ 1 ) 

1 2T'+1 / -■ \T' l 2W+l ( ^\W+T'-T 
Y 2J"+1\ L ) y 2T+1 ^ ' 

v (_i\v- t i r <J"M" r TM T r J'M' 

/ 2T'+1 ^ -i\T' / 2W+1 / i\W+r'-T 

Y2J 77 TT^ _i ^ Y^+ r ^~ ' 

x(-lf+^''+''V2J^V2J^^ ^ £ ](133) 
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see Eq. 8.7(15) in Ref. [6|. This gives: 

p:% u : w ' J " M \r) = {-if £ 6 3+ ^v^<§iki|i> 



Nlj 
N'l'j' 



x E ^^Wz(^)^Tr(^illP J ''ll0iv'ZY) 



k' 



X 



( _ 1)J '+i^+i-72 v ^JTTv / 2jTT v / 2j 7 TT 
t> I V T 



/ (2fc+l)(2fc'+l) r T0 y- (n ,X 
/ / V 1^(27+1) X ^M^'O^Mt^, </>J 

x (-i) J + fe - t/, - r v /(2A; + l)(2fc' + 1)(2T' + 1)(2T + 1) 
f I k U ) 
x\ V k' U> 
[ V T W J 

/ 2T'+1 I i\T' /2W±1{ i\W+T'-Tf i\T'+W+J"+J' 

Finally we have: 

~uUu'U'W,J" M" i \ , n i +t ,, 3 +«+«' -(6r) 2 /l|, Mix IT 

PvJ'M' {') — \ l ) e \2ll a «ll2/ Y 4tt 

JVijfc 
N'l'j' k' 



T'TMt 

' (2T'+l)(2T'+l)(2k+l)(2k'+l)(2J'+l)(2W+l)(2j+l)(2j'+l) 
X v (2T+1) 



j y </" ] r i k u 

11 u > < V k! V 



v V J" 



22 \\ (\TTW 

I V V T T W 
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X C T J^M'J" M"YtMt(®> 0) j 



(135) 



where we have also used Eq. (|9~Uj) . This gives expression (I98p and definitions 
of radial form factors (j9"9"|) and coefficients fllOOp . 

C. Derivation of Eq. ( 1X131 ) 

After inserting Eq. (11 111) into Eq. f 1 1 1 2 j) we must sum up products of four 
Clebsh-Gordan coefficients, similarly as in Eq. (11291) . that is, 



E r 2 m 's r jm 3 ^j'mr j'm'- 
1 , , , 1 °im,/"M" ( - y ,, ,1 , 



= (-l) M "~ m 'L- l ~ v 'y/2^2f + ly/2f + l^/2]Tl 

( J f I" ] 

X ^Vm'l-mPv'-l'V'M 1 ' \ 2 2 V ( ' ( 136 ) 

( / I' T'J 

and then as in Eq. f 1 1 3 1 j) : 

El _i\-rn', r km k r k'm' k n T'M' T r k'm' k 

\ L ) ^lm l IM I ^Vrn[I'-M' I ( -'Vrn' l l-rnC-'TM T krn k 

mim[m k m' k 



x £ Hr^tA ( ' * M ■ (137) 



wm^ \ T T W 



We can now perform the summation over and Mj, which gives the 
factor (— 1) /+/ ~ L Sliw'Sm' m' an d allows for a summation over W and M(^. 
After inserting all these results into Eq. f 1 1 X 2 [) . we obtain 

(^ii^''(p j '' m '')ii^)=E E E E 

M" n'L'v'J' mlm'P M' L p! 

E j, M < j-iy'-M' ( U m' K n'L' i -, 



, ^ L' M' L v' fi' y/2J'+l ^ ± ^mlm'l'\ i 

M' 



s/ \^\^V^ 1 l|l\ / (2T+l)(2fc+l) /-yfc'O 

TMt k> k 
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x 

X 



X 



Lj J',-M'J"M" 

(-l)^-'V2V27TIV2j+l 
f J f I" 

A '«' -li'V'M" \ 2 2 y 

t'm' t [ i v r 

(-l) M i(-l) i + fc + fc '- / ' v /(2T' + l) 
Z /v / 

*' # ^ x-i) 7+y '^'(i38) 

T' T L' 



Now, we have to sum up products of three Clebsh-Gordan coefficients, simi- 
larly as in Eq. (11331) . that is, 



\ L ) Ly L'M' L v' l i' Ly v'-tM'I"M' I '^'[ 
M' L M' T fi' 

_ I ZT'+K i\T' l 2L'+l ( i \L'+T'-T+M T 

This gives: 
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Finally, summations over M' and My give the factor Sj"j"5m'/m"j that is: 
(4>N'vA\W(? ,ua )\\<l>mj) 

= 5 I „j„(-l) J ^2 K mlL'I'{- l Y (\\\ a A\\)\[i 

n'L'v'J' mlra'V 

£(-lf (2T + 1) f rMr6 3+ ™'Fr, /w ,,(6r)^i;, J ,(6r)F^(6r)e- 2 ^ 2 



x 

Tkk' 



Y / -l\k( Y\T( i\I-L' r T0 / (2fc+l)(2fc'+l)(2L'+l)(2j+l)(2j'+l) 
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where we have also used Eq. (1901) . This shows explicitly that in the field 
calculated for the density matrix with multipolarity J", only the multipole 
J" appears. 
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